module projectors
  !
  ! Projection Generator
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
#define green(i) (iand(mana,i)/=0)
  use ground_control
  use oct_pet
  implicit none

  private
  public initProjectors
  public bproj
  public fproj
  public enableSubsets,disableSubsets
  public toNextSubset
  public setCurrentSubset

  !
  integer nmaxth

  !
  real,target,allocatable::x_js(:,:,:) ! backprojected image
  real,pointer::xj_in(:,:)
  real,pointer::xj_out(:,:)
  real,pointer::pi_in(:,:,:,:,:,:)
  real,pointer::pi_out(:,:,:,:,:,:)
  integer num_of_pi
  integer,allocatable::thmap(:)
  integer num_of_xjs

  !
  integer,allocatable,target::mapToRjs(:,:)

  type t_projargs
     type(t_sysmat),pointer::sm
     integer,pointer::it1,it2,da
     integer pid,gid
     integer npi_4
  end type t_projargs
  type(t_projargs),allocatable,target::pas(:)

  integer,allocatable::dt1s(:)
  integer,allocatable::da1s(:)  
  integer,public::nos=1
  integer,public::mos=1
  integer,public::current_sid

  external setma
  pointer(psetma,setma)
  
  external sump
  pointer(psump,sump)

  type t_ordered_subset
     integer n
     integer,allocatable::idt(:)
  end type t_ordered_subset
  type(t_ordered_subset),allocatable::ost(:)

contains

  subroutine echo(mess)
    use suty
    character*(*),intent(in)::mess
    character(MAXPROCLEN),save::theProc
    character(MAXPROCLEN) proc
    if(.not.green(OPT_VERBOSE)) return
    proc=getProcname(mess)
    if(proc=="") then
       proc=theProc
    else
       theProc=proc
       proc=" "
    end if
    call echo_green(mess,"","pr")
  end subroutine echo

  !
  subroutine initProjectors(penv)
    use omp_lib
    use suty
    type(t_opet_env),intent(in)::penv
    integer j,jr,jjt,jt,rot

    nmaxth=omp_get_max_threads()
    if(allocated(thmap)) deallocate(thmap)
    allocate(thmap(nmaxth))

    call setSubsetsOrder
    call echo("initProjectors: number of subsets: "//trim(i2a(nos)))

    if(allocated(mapToRjs)) deallocate(mapToRjs)
    allocate(mapToRjs(nj,0:md-1))
    j=1
    mapToRjs(1,:)=1
    do jr=1,nr
       do jt=1,maxjts(jr)
          j=j+1
          do rot=0,md-1
             jjt=jt+rot*jr*ndt
             if(jjt>maxjts(jr)) jjt=jjt-maxjts(jr)
             mapToRjs(j,rot)=jrt2j(jr,jjt)
          end do
       end do
    end do

    if(allocated(x_js)) deallocate(x_js)
    allocate(x_js(nj,nz,nmaxth))
    num_of_xjs=nj*nz*nmaxth

    if(allocated(pas)) deallocate(pas)
    allocate(pas(nmaxth))

    num_of_pi=ma*mt*ma*mt*npid*ngid

  contains

    subroutine setSubsetsOrder
      integer is

      if(allocated(dt1s)) deallocate(dt1s)
      if(allocated(da1s)) deallocate(da1s)      
      current_sid=1
      nos=1
      psump = loc(sumc)      
      
      if(penv%number_of_subsets>1) then
         mos=penv%number_of_subsets
         nos=penv%number_of_subsets * 2
      else
         return
      end if
      psump = loc(sums)
      allocate(da1s(nos))
      allocate(dt1s(nos))
      do is=1,mos
         da1s(is)=is-1
         dt1s(is)=-maxcdt+(is-1)
      end do
      do is=mos+1,nos
         da1s(is)=da1s(is-mos)
         dt1s(is)=dt1s(nos-is+1)
      end do
    end subroutine setSubsetsOrder

  end subroutine initProjectors

  !
  subroutine fproj(x_j,p_i)
    use suty
    real,intent(in),target::x_j(:,:)
    real,intent(out),target::p_i(:,:,:,:,:,:)    

    ! 
    call clemem8(loc(p_i),num_of_pi)
    xj_in => x_j
    pi_out => p_i
    psetma = loc(setpi)
    call runsum

  end subroutine fproj

  !
  subroutine bproj(p_i,x_j)
    use suty
    real,intent(in),target::p_i(:,:,:,:,:,:)
    real,intent(out),target::x_j(:,:)
    integer i

    !
    call clemem8(loc(x_js),num_of_xjs)
    thmap=0
    pi_in => p_i
    psetma = loc(setdj)
    call runsum
    
    do i=2,nmaxth
       if(thmap(i)==0) cycle       
       x_js(:,:,1)=x_js(:,:,1)+x_js(:,:,i)
    end do
    x_j=x_js(:,:,1)

  end subroutine bproj

  subroutine runsum
    integer i
    !$omp parallel do
    do i=1,ngid*npid
       call sumc(i)
    end do
    !$omp end parallel do    
  end subroutine runsum
  
  subroutine sumc(pgid)
    use omp_lib
    integer,intent(in)::pgid
    integer ith
    integer,target::it1,it2,da,dt
    integer id1,id2,idd
    type(t_projargs),pointer::pa
    type(t_block_pair)::idts(*)
    pointer(pidts,idts)
    ith=omp_get_thread_num()+1      
    thmap(ith)=1
    pa => pas(ith)
    pa%da  => da
    pa%it2 => it2
    pa%it1 => it1

    pidts = loc(idt)
    id1=abs(idts(pgid)%id1)
    id2=idts(pgid)%id2
    idd=abs(4-(id2-id1))+1    

    pa%npi_4=idts(pgid)%rot
    pa%pid=abs(ridt(id1,id2)%pid)
    pa%gid=abs(ridt(id1,id2)%gid)

    do da=0,maxcda
       do dt=-maxcdt,maxcdt
          do it1=max(1-dt,1),min(mt-dt,mt)
             it2=dt+it1
             pa%sm => sysmat(it1,it2,da+1,idd)             
             if(pa%sm%n==0) cycle
             call setma(ith)
          end do
       end do
    end do
  end subroutine sumc

  subroutine setpi(i)
    integer,intent(in)::i
    type(t_projargs),pointer::pa
    type(t_sysmat),pointer::sm
    integer ia1
    integer jz,j,djz
    integer ia1r,ia2r,ia2
    integer,pointer::Rjs(:)
    integer*4 js(*),jzs(*)
    real,pointer::theProj,theProjR
    real*8 ws(*)
    pointer(pj,js)
    pointer(pjz,jzs)
    pointer(pw,ws)
    integer k

    pa => pas(i)
    sm => pa%sm

    pj=sm%p
    pjz=pj+sm%n*4
    pw=pjz+sm%n*4

    Rjs => mapToRjs(:,pa%npi_4)

    do ia1=1,ma-pa%da

       djz=(ia1-1)*nint(cda/dz)

       ia2=ia1+pa%da
       ia1r=(ma-ia1)+1
       ia2r=(ma-ia2)+1
              
       theProj  => pi_out(pa%it1,ia1, pa%it2,ia2, pa%pid,pa%gid)
       theProjR => pi_out(pa%it1,ia1r,pa%it2,ia2r,pa%pid,pa%gid)
       
       do k=1,sm%n
          jz=jzs(k)+djz
          j=Rjs(js(k))
          if(xj_in(j,jz)==0) cycle
          theProj=theProj+xj_in(j,jz)*ws(k)
       end do

       if(pa%da==0) cycle
       do k=1,sm%n
          jz=nz-(jzs(k)+djz)+1
          j=Rjs(js(k))
          if(xj_in(j,jz)==0) cycle
          theProjR=theProjR+xj_in(j,jz)*ws(k)
       end do
       
    end do

  end subroutine setpi

  subroutine setdj(i)
    integer,intent(in)::i
    type(t_projargs),pointer::pa
    type(t_sysmat),pointer::sm
    integer ia1
    integer jz,j,djz
    integer ia1r,ia2r,ia2
    real,pointer::theProj,theProjR
    real,pointer::x_j(:,:)
    integer,pointer::Rjs(:)
    integer*4 js(*),jzs(*)
    real*8 ws(*)
    pointer(pj,js)
    pointer(pjz,jzs)
    pointer(pw,ws)
    integer k

    pa => pas(i)
    sm => pa%sm

    pj=sm%p
    pjz=pj+sm%n*4
    pw=pjz+sm%n*4

    Rjs => mapToRjs(:,pa%npi_4)

    x_j => x_js(:,:,i)

    do ia1=1,ma-pa%da

       djz=(ia1-1)*nint(cda/dz)

       ia2=ia1+pa%da
       ia1r=(ma-ia1)+1
       ia2r=(ma-ia2)+1      
       
       theProj  => pi_in(pa%it1,ia1, pa%it2,ia2, pa%pid,pa%gid)
       theProjR => pi_in(pa%it1,ia1r,pa%it2,ia2r,pa%pid,pa%gid)

       if(theProj/=0) then
          do k=1,sm%n
             jz=jzs(k)+djz
             j=Rjs(js(k))
             x_j(j,jz)=x_j(j,jz)+theProj*ws(k)
          end do
       end if

       if(pa%da==0.or.theProjR==0) cycle
       do k=1,sm%n
          jz=(nz-(jzs(k)+djz))+1
          j=Rjs(js(k))
          x_j(j,jz)=x_j(j,jz)+theProjR*ws(k)
       end do

    end do
  end subroutine setdj

!!!!
  
  subroutine setCurrentSubset(is)
    integer,intent(in)::is
    current_sid=is
  end subroutine setCurrentSubset

  subroutine toNextSubset
    current_sid = current_sid + 1
    if(current_sid > nos) current_sid = 1
  end subroutine toNextSubset
  
  subroutine disableSubsets
    call enableSubsets(.false.)
  end subroutine disableSubsets
  
  subroutine enableSubsets(on)
    logical,intent(in),optional::on
    logical on_
    if(present(on)) then
       on_=on
    else
       on_=.true.
    end if
    on_ = on_ .and. (nos>1)
    if(on_) then
       psump = loc(sums)
       current_sid = 1
    else
       psump = loc(sumc)       
    end if
  end subroutine enableSubsets

  subroutine sums(pgid)
    use omp_lib
    integer,intent(in)::pgid
    integer ith
    integer,target::it1,it2,da,dt
    integer id1,id2,idd
    type(t_projargs),pointer::pa
    type(t_block_pair)::idts(*)
    pointer(pidts,idts)

    ith=omp_get_thread_num()+1      
    thmap(ith)=1
    pa => pas(ith)
    pa%da  => da
    pa%it2 => it2
    pa%it1 => it1

    pidts = loc(idt)
    id1=abs(idts(pgid)%id1)
    id2=idts(pgid)%id2
    idd=abs(4-(id2-id1))+1    

    pa%npi_4=idts(pgid)%rot
    pa%pid=abs(ridt(id1,id2)%pid)
    pa%gid=abs(ridt(id1,id2)%gid)

    do da=da1s(current_sid),maxcda,mos
       do dt=dt1s(current_sid),maxcdt,mos
          do it1=max(1-dt,1),min(mt-dt,mt)
             it2=mt-(dt+it1)+1
             pa%sm => sysmat(it1,it2,da+1,idd)             
             if(pa%sm%n==0) cycle
             call setma(ith)
          end do
       end do
    end do

  end subroutine sums

end module projectors
